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We have studied the d-wave pairing-instability in the two-dimensional Holstein-Hubbard model 
at the level of a full fluctuation exchange approximation which treats both Coulomb and electron- 
phonon (EP) interaction diagrammatically on an equal footing. A generalized numerical renormal- 
ization group technique has been developed to solve the resulting self-consistent held equations. The 
d-wave superconducting phase diagram shows an optimal T c at electron concentration (n) ~ 0.9 for 
the purely electronic Hubbard system. The EP interaction suppresses the d-wave T c which drops to 
zero when the phonon-mediated on-site attraction U p becomes comparable to the on-site Coulomb 
repulsion U. The isotope exponent a is negative in this model and small compared to the classi- 
cal BCS value ubcs = | or compared to typical observed values in non-optimally doped cuprate 
superconductors. 



In recent years, growing experimental evidence has 
suggested that YBa2Cu3C>7 and, possibly, other cuprates 



are d^.^-superconductorsEJ. Anti-ferromagnetic (AF) 
spin fluctuation (SF) exchange has been proposed as a 
possible candidate mechanism for d-wave pairing. These 
AF spin fluctuation models are based on the notion that 
short-range, dynamical AF spin correlations, caused by 
the strong local Coulomb repulsion in the cuprites, may 
lead to a spatially extended pairing attraction.Erta Start- 
ing from purely electronic models, such as the Hubbard 
Hamiltonian, coupling to lattice vibrational degrees of 
freedom is usually neglected in this picture. However, ex- 
cept near certain "optimal" doping concentrations, many 
cuprates, including YBa2Cu307, exhibit a quite notice- 
able doping dependent isotope effecttj. This indicates 
that electron-phonon (EP) interactions could be impor- 
tant and should be included in the theory. 

The goal of the present paper is to study the com- 
petition between phonons and AF spin fluctuation ex- 
change by means of a self-consistent diagrammatic ap- 
proach which explicitly includes phonon renormalizations 
to the AF spin fluctuations at the level of the effective 
interaction vertices. We formulate a full fluctuation ex- 
change (FLEX)Elu approximation which treats Coulomb 
and EP contributions to the electron-electron interac- 
tion potential entirely on an equal footing. Qur work 
goes substantially beyond previous treatmentsQu which 
have included phonon effects only at the level of the one- 
particle self-energy. 

Due to the retarded nature of the EP interaction, the 
problem is numerically not directbz amenable to the re- 
cent fast Fourier transfomi (FFTp or numerical renor- 
malization group (NRGp methods, developed for for 
the FLEX approximation to the pure Hubbard model. 



The present FLEX equations require certain large-scale 
fcrmion frequency matrix inversions which are numeri- 
cally about 4 orders of magnitude more demanding than 
FLEX calculations with short-range instantaneous inter- 
actions. The numerical solution of this problem can be 
achieved only by means of a generalized, highly efficient 
matrix version of the original NRG methods which we 
have developed. 

We start from the simplest microscopic Hamiltonian 
which includes both an on-site Hubbard U Coulomb re- 
pulsion and a local EP coupling to an, Einstein phonon 
branch, the Holstein-Hubbard modell 
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with a nearest neighbor hopping t, chemical potential 
/i, on-site Coulomb repulsion U , on-site EP coupling 
constant C, force constant K, and ionic oscillator mass 
M. The ct (c i(7 ) is the electron creation (annihila- 
tion) operator at site i and spin a; is the number 
operator; and is the local ionic displacement at lat- 
tice site i. The dispersionless bare phonon frequency is 
Qq = (K/M) 1 / 2 and the phonon-mediated on-site at- 
traction is Up = C 2 /K. n 

Previous self-consistent field (SCF) studieso of the 
Holstein-Hubbard system have ignored the electron- 
electron exchange scattering which arises from the Pauli 
exclusion principle. The importance of this exchange 
vertexoQ can be most easily demonstrated in the limit 
of the negative-U Hubbard model with U — —\U\ and 
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Up = 0. In this case, the direct interaction will give 
2U (after summing over the electron spin index) and 
the exchange interaction contributes —U. The simplest 
mean field theory will then predict the CDW instability 
to occur at \U\x P h(T) = 1 (with exchange interaction) 
instead of 2\U\x P h(T) = 1 (without exchange interac- 
tion). In the positive-!/ Hubbard model, the exchange 
interaction enhances the spin fluctuations and thus helps 
the d-wave instability while at the same time weaken- 
ing the charge fluctuations. Also, high phonon frequen- 
cies or a flat electron band near the Fermi surface will 
tend to enhance the effect of the exchange vertex. To 
study particle-hole and particle-particle instabilities in 
this model, it is thus necessary to include both Coulomb 
and phonon contributions to the exchange vertex. 

The bare interaction vertices for the FLEX equations 
shown in Fig.jl], include the particle-hole [Fig. |](b) and 
(c)] and the particle-particle [Fig. 0(d)] bare vertices, 
due to both the Hubbard U and the phonon propagator 
v p (iu m ) = -U p nl/{nl + i^) [Fig. 0(a)] for boson Mat- 
subara frequeruy- 1 i' TO = 2mnT. The one-particle self- 
energy is then:! 
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R n , ni (Ak) = < ni (Ao;) x 



fori?= ^° r f ' (6) 
1 X P p{Ak;iw ni ) 1 b or T , 

where k = (k, iu) n ), Ak = (k — k', iw n — iu) n '), Lj n 
2 a + V)ttT, the Green's function G(k) — [iu> n — £k — 
£&] , and the tight binding band = — 2i(cos k x + 
cos k y ) — \i. The v^ ni are the bare vertices shown in 
Fig. l(b-d). The bare particle-hole and particle-particle 
fluctuation functions are defined as: 
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FIG. 1. The bare vertices of the Holstein-Hubbard model 
in the FLEX approximation, (a) The interactions of the 
on-site Coulomb U and Einstein phonon v p (iv m ). (b) Den- 
sity vertex v° x ^ ni (iv m ) = [2v p (iu m ) — v p (iuj ni — iw nA ) + U] 



(c) Magnetic vertex v„ l!ni ( 
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X ph (q; iuj n ) = ~^Y1 G ^ k + (? ) G ( fc ) ' 
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X PP (q;iu n )= -^G(fc + <z)G(-fc) . (8) 



Because of the retarded nature of v p (w m ) 1 the bare ver- 
tices in Fig.l(b-d) depend explicitly on the internal fre- 
quency transfer and iuj n can not be summed out here. 
The numerically most challenging part of the SCF calcu- 
lation is thus the evaluation of the fluctuation potentials, 
V ph and V pp , because of the required fermion frequency 
matrix inversion in Eqs. (4,5). In a brute-force approach, 
this matrix dimension can become as large as 500 x 500 
(the size of the entire fermion Matsubara frequency set) 
at T c . Recent FFT and NRG techniques for the pure 
Hubbard FLEX equations are not directly applicable and 
we had to develop a generalized " fermion_matrix" NRG 
method to handle the numerics efficiently^. n 

Building upon the original NRG methocH, our gen- 
eralized fermion matrix NRG employs a frequency-RG 
operation which separates the frequency space into high 
and low regions at each temperature in such a way that 
fermion frequency matrix dimensions will not increase as 
fast as T _1 when T is lowered. For a typical 8x8 ma- 
trix dimension, starting at some high temperature and 
large fermion frequency cut-off, we can keep the dimen- 
sion below 30 x 30 after 6 factor-2 frequency-RG steps, 
corresponding to a temperature reduction from 4.0i to 
0.0625£. The fermion cut-off in the frequency RG steps 
changes from lOOi (~ 12x bandwidth) to 1.6* (~ 50xT c ). 
We then keep the same cut-off (without any further RG 
operation) and decrease the temperature slowly until the 
d-wave instability occuraj. A typical matrix dimension 
is then about 50 x 50 near the d-wave instability. The 
d-wave instability is determined by solving the eigen- 
value problem for the singlet particle-particle kernel, con- 
structed from the full interaction potential in Fig. |l|(a), 
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FIG. 2. Phase diagram of the d-wave instability with 
U/t — 4 and Qo/t — 0.5 and different U p on a 32 2 lattice. 
Solid line represents a Hubbard model calculation using the 
NRG method in Ref. [9] which does not involve the matrix 
inversion in Eqs. (4,5). 

following Ref. |^. T c is reached when the maximum paring 
eigenvalue Kd(T) = 1. A 32x32 k-mesh covering the 
full 1st Brillouin zone has been employed in all of the 
calculations reported below. 

Figure || shows the ci-wave superconducting phase dia- 
gram of the Holstein-Hubbard model with intermediate 
Hubbard U/t = 4 for various EP coupling U p = 0, 2, 
and 3 and an Einstein phonon frequency £lo/t = 0.5. 
Increasing the electron concentration (n) towards half- 
filling ((n) ~ 1) initially enhances the ci-wave T c un- 
til it reaches a maximum around (n) ~ 0.9. Beyond 
that point, at hole doping x = 1 — (n) below ~ 8%, 
the strong AF fluctuations actually reduce the ci-wave 
T c . This should be contrasted with early Hubbard model 
FLEX results in Ref. |[ where the detailed shape of the 
magnetic-superconducting boundary was not well deter- 
mined and the ci-wave T c calculation was stopped when 
the magnetic eigenvalue exceeded unity. Here, we have 
used a finer lattice mesh and larger cut-off frequency and 
carefully pushed the ci-wave calculations toward smaller 
hole doping. 

In order to get a deeper understanding of the ori- 
gin for the T c maximum in pEig. 1, we have carried 
out a McMillan-type analysisEJ of the underlying pair- 
ing equationsel by estimating the dimensionless Eliash- 
berg parameter Ad, which measures the pairing poten- 
tial strength averaged over the Fermi surface in the 
ci-wave channel, and Xz = — 9 w RcE(fc)| aJ= j + , which 
measures the strength of the quasi-particlc mass en- 
hancement, as well as the pair-breaking strength 7 = 

ImE(fc) /T c , due to the quasi-particle damping. 

uj=i0+ I — I 

Expressed in terms of " renormalized" parameters^ A d = 
A<j/ (1 + Az) and 7* = 7/(1 + Xz), our results show that 



FIG. 3. d-wave T c vs. EP potential U p for U/t = 4, 
(n) = 0.9, and flo/t = 0.5. Inset is the maximum d-wave 
eigenvalue vs. T for U p /t — 0, 1,2,3, and 4 (from top to bot- 
tom). T c for Up/t = 4 is extrapolated from the Hid data at the 
lowest available T. 

both the pairing strength A^ and the pair-breaking 
strength 7* are monotonically increasing as (n) is pushed 
towards half-filling. Thus there are (at least) two com- 
peting effects at work here: While raising AJ increases 
T c , raising 7* decreases it. Apparently, for overdoping, 
the doping variation of the pairing strength A^ domi- 
nates T c , making T c initially rise with increasing (n). On 
the other hand, for underdoping, close to half-filling, the 
doping variation of the pair-breaking strength 7* domi- 
nates and causes T c to decrease with increasing (n). An 
additional, related effect is that the overall AF spin fluc- 
tuation energy scale softens as (n) approaches half-filling. 
This lowering of the relevant "boson" energy scale will 
also lower T c . 

The primary effect of EP coupling is to suppress the ci- 
wave T c , shown in Fig. || as a function of the EP potential 
strength U p (from to 4t) at fixed U/t = 4, electron fill- 
ing (n) = 0.9, and Einstein phonon frequency £lo/t = 0.5. 
The corresponding maximum ci-wave pairing eigenvalues 
Kd as a function of temperature are plotted in the in- 
set. Note that the ci-wave T c drops to "almost zero" (i.e. 
numerically inaccessible values) when U p becomes com- 
parable to r£/. This behavior is different from our earlier 
calculationQ which ignored the phonon renormalization 
of the bare interaction vertices. In that cased, T c was 
suppressed only by the EP self-energy contribution, the 
suppression was much more gradual and T c dropped only 
by about one half between U p = to U p ~ U. Here, by 
contrast, the EP interaction directly counteracts the on- 
site Coulomb repulsion and thereby suppresses the AF 
spin fluctuation mediated pairing potential. The ci-wave 
Eliashberg pairing strength A^ and the pair-breaking 
strength 7 are indeed found to be strongly suppressed 
by the EP interaction. 
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TABLE I. Isotope exponent a for U/t = 4 and U p /t = 2. 
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An important feature of EP coupling is that it in- 
troduces an isotope effect into the electronic d-wave 
pairing mechanism. Table I shows results for the iso- 
tope exponent a — — d log T c /d log M which becomes 
a = id log T c /d log f2o in the present model, since the 
isotopic mass M enters only through Qq. In our previ- 
ous studieso of the isotope effect, where the EP effect 
on the magnetic bare vertices was neglected, the iso- 
tope exponent a was quite small and negative for re- 
alistic phonon energies f2o- Here, as shown in Table I, we 
find qualitatively the same result, even though the EP 
effect on the AF spin fluctuations is now explicitly taken 
into account and suppresses T c much more strongly. It 
is interesting to note that the absolute value of a has 
a minimum at the optimal doping concentration, a fea- 
ture qualitatively reminiscent of the doping dependent 
isotope data in many cuprate systemsla However the ob- 
served overall magnitude of the effect in non-optimally 
doped cuprates, |a| ~ 0.5 — lja is much larger than the 
present model predicts. This finding further supports the 
notion that the EP coupling in the cuprates could be ef- 
fectively very much enhanced cpmpared to conventional 
"strong-coupling" EP systems. 

In conclusion, we have studied the d-wave supercon- 
ducting instability of the Holstein-Hubbard model in the 
FLEX approximation by means of a generalized, ma- 
trix version of the numerical renormalization group tech- 
nique. Upon including both the particle-particle and 
particle-hole fluctuations, the d-wave T c shows a max- 
imum at electron filling (n) ~ 0.9. The d-wave T c is 
suppressed by increasing the EP potential U p and tends 
to zero when U p ~ U. Finally, the isotope exponents 
a are negative and small in magnitude, with typically 
\a\ < 0.1. While \a\ exhibits a minimum at optimal dop- 
ing concentration, the overall magnitude is far too small 
to explain observed isotope data in the cuprates. Our full 
FLEX results support the conclusions of earlier isotope 
calculations by the present authors .0 
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